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Abstract 



In this article, the extensional flow and viscosity and the converging-diverging 
geometry were examined as the basis of the peculiar viscoelastic behavior in porous 
media. The modified Bautista-Manero model, which successfully describes shear- 
thinning, elasticity and thixotropic time-dependency, was used for modeling the 
flow of viscoelastic materials which also show thixotropic attributes. An algorithm, 
originally proposed by Philippe Tardy, that employs this model to simulate steady- 
state time-dependent flow was implemented in a non-Newtonian flow simulation 
code using pore-scale modeling and the initial results were analyzed. The findings 
are encouraging for further future development. 
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1 Introduction 

Non-Newtonian fluids are commonly divided into three broad groups [1, 2]: 

1. Time-independent: those fluids for which the strain rate at a given point is 
solely dependent upon the instantaneous stress at that point. 

2. Viscoelastic: those fluids that show partial elastic recovery upon the removal 
of a deforming stress. 

3. Time-dependent: those fluids for which the strain rate is a function of both 
the magnitude and the duration of stress and possibly of the time lapse 
between consecutive applications of stress. 

A large number of models have been proposed in the literature to model all 
types of non-Newtonian fluids under various flow conditions. Most these models 
are basically empirical in nature and arising from curve-fitting exercises [3]. In 
this article, we investigate a non-Newtonian rheological model (namely Bautista- 
Manero model) that can be used to describe viscoelastic behavior with thixotropic 
attributes in the context of pore-scale modeling of non-Newtonian flow through 
porous media. 
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2 Viscoelastic Fluids 

Viscoelastic substances exhibit a dual nature of behavior by showing signs of both 
viscous fluids and elastic solids. In its most simple form, viscoelasticity can be 
modeled by combining Newton's law for viscous fluids (stress oc rate of strain) 
with Hook's law for elastic solids (stress oc strain), as given by the original Maxwell 
model and extended by the Convected Maxwell models for the nonlinear viscoelastic 
fluids. Although this idealization predicts several basic viscoelastic phenomena, it 
does so only qualitatively [4]. 

Polymeric fluids often show strong viscoelastic effects, which can include shear- 
thinning, extension thickening, viscoelastic normal stresses, and time-dependent 
rheology phenomena. The equations describing the flow of viscoelastic fluids consist 
of the basic laws of continuum mechanics and the rheological equation of state, 
or constitutive equation, describing a particular fluid and relates the viscoelastic 
stress to the deformation history. The quest is to derive a model that is as simple 
as possible, involving the minimum number of variables and parameters, and yet 
having the capability to predict the viscoelastic behavior in complex flows [5] . 

No theory is yet available that can adequately describe all of the observed vis- 
coelastic phenomena in a variety of flows. However, many differential and integral 
viscoelastic constitutive models have been proposed in the literature. What is com- 
mon to all these is the presence of at least one characteristic time parameter to 
account for the fluid memory, that is the stress at the present time depends upon 
the strain or rate-of-strain for all past times, but with an exponentially fading 
memory [6, 7, 8, 9, 10]. 

The behavior of viscoelastic fluids is drastically different from that of Newtonian 
and inelastic non-Newtonian fluids. This includes the presence of normal stresses 
in shear flows, sensitivity to deformation type, and memory effects such as stress 
relaxation and time-dependent viscosity. These features underlie the observed pe- 
culiar viscoelastic phenomena such as rod-climbing (Weissenberg effect), die swell 
and open-channel siphon [4, 11]. Most viscoelastic fluids exhibit shear-thinning and 
an elongational viscosity that is both strain and extensional strain rate dependent, 
in contrast to Newtonian fluids where the elongational viscosity is constant [11]. 

The behavior of viscoelastic fluids at any time is dependent on their recent 
deformation history, that is they possess a fading memory of their past. Indeed 
a material that has no memory cannot be elastic, since it has no way of remem- 
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bering its original shape. Consequently, an ideal viscoelastic fluid should behave 
as an elastic solid in sufficiently rapid deformations and as a Newtonian liquid in 
sufficiently slow deformations. The justification is that the larger the strain rate, 
the more strain is imposed on the sample within the memory span of the fluid 
[4, 11, 12]. 

Many materials are viscoelastic but at different time scales that may not be 
reached. Dependent on the time scale of the flow, viscoelastic materials mainly 
show viscous or elastic behavior. The particular response of a sample in a given 
experiment depends on the time scale of the experiment in relation to a natural 
time of the material. Thus, if the experiment is relatively slow, the sample will 
appear to be viscous rather than elastic, whereas, if the experiment is relatively 
fast, it will appear to be elastic rather than viscous. At intermediate time scales 
mixed viscoelastic response is observed. Therefore the concept of a natural time 
of a material is important in characterizing the material as viscous or elastic. The 
ratio between the material time scale and the time scale of the flow is indicated by 
a non-dimensional number: the Deborah or the Weissenberg number [3, 13]. 

A common feature of viscoelastic fluids is stress relaxation after a sudden shear- 
ing displacement where stress overshoots to a maximum then starts decreasing 
exponentially and eventually settles to a steady state value. This phenomenon 
also takes place on cessation of steady shear flow where stress decays over a finite 
measurable length of time. This reveals that viscoelastic fluids are able to store 
and release energy in contrast to inelastic fluids which react instantaneously to the 
imposed deformation [4, 10, 12]. 

A defining characteristic of viscoelastic materials associated with stress relax- 
ation is the relaxation time which may be defined as the time required for the 
shear stress in a simple shear flow to return to zero under constant strain condi- 
tion. Hence for a Hookean elastic solid the relaxation time is infinite, while for a 
Newtonian fluid the relaxation of the stress is immediate and the relaxation time 
is zero. Relaxation times which are infinite or zero are never realized in reality 
as they correspond to the mathematical idealization of Hookean elastic solids and 
Newtonian liquids. In practice, stress relaxation after the imposition of constant 
strain condition takes place over some finite non-zero time interval [8]. 

The complexity of viscoelasticity is phenomenal and the subject is notorious 
for being extremely difficult and challenging. The constitutive equations for vis- 
coelastic fluids are much too complex to be treated in a general manner. Further 
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complications arise from the confusion created by the presence of other phenomena 
such as wall effects and polymer-wall interactions, and these appear to be system 
specific. Therefore, it is doubtful that a general fluid model capable of predicting 
all the flow responses of viscoelastic system with enough mathematical simplicity 
or tractability can emerge in the foreseeable future [10, 14, 15]. Understandably, 
despite the huge amount of literature composed in the last few decades on this 
subject, almost all these studies have investigated very simple cases in which sub- 
stantial simplifications have been made using basic viscoelastic models. 

In the last few decades, a general consensus has emerged that in the flow of 
viscoelastic fluids through porous media elastic effects should arise, though their 
precise nature is unknown or controversial. In porous media, viscoelastic effects 
can be important in certain cases. When they are, the actual pressure gradient 
will exceed the purely viscous gradient beyond a critical flow rate, as observed 
by several investigators. The normal stresses of high molecular polymer solutions 
can explain in part the high flow resistances encountered during viscoelastic flow 
through porous media. It is argued that the very high normal stress differences and 
Trouton ratios (defined as extensional to shear viscosity) associated with polymeric 
fluids will produce increasing values of apparent viscosity when flow channels in 
the porous medium are of rapidly changing cross section. 

Important aspects of non-Newtonian flow in general and viscoelastic flow in 
particular through porous media are still presenting serious challenge for modeling 
and quantification. There are intrinsic difficulties of characterizing non-Newtonian 
effects in the flow of polymer solutions and the complexities of the local geometry 
of the porous medium which give rise to a complex and pore-space-dependent 
flow field in which shear and extension coexist in various proportions that cannot 
be quantified. Flows through porous media cannot be classified as pure shear 
flows as the converging- diverging passages impose a predominantly extensional flow 
fields especially at high flow rates. Moreover, the extension viscosity of many non- 
Newtonian fluids increases dramatically with the extension rate. As a consequence, 
the relationship between the pressure drop and flow rate very often do not follow 
the observed Newtonian and inelastic non-Newtonian trend. Further complication 
arises from the fact that for complex fluids the stress depends not only on whether 
the flow is a shearing, extensional, or mixed type, but also on the whole history of 
the velocity gradient [5, 16, 17, 18, 19, 20]. 
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3 Important Aspects for Viscoelastic Flow in Porous 
Media 

Strong experimental evidence indicates that the flow of viscoelastic fluids through 
packed beds can exhibit rapid increases in the pressure drop, or an increase in 
the apparent viscosity, above that expected for a comparable purely viscous fluid. 
This increase has been attributed to the extensional nature of the flow field in the 
pores caused by the successive expansions and contractions that a fluid element 
experiences as it traverses the pore space. Even though the flow field at pore level is 
not an ideal extensional flow due to the presence of shear and rotation, the increase 
in flow resistance is referred to as an extension thickening effect [19, 21, 22, 23]. 

There are two major interrelated aspects that have strong impact on the flow 
through porous media. These are extensional flow and converging-diverging geom- 
etry. 

3.1 Extensional Flow 

One complexity in the flow in general and through porous media in particular 
usually arises from the coexistence of shear and extensional components; sometimes 
with the added complication of inertia. Pure shear or elongational flow is very much 
the exception in practical situations, especially in the flow through porous media. 
By far the most common situation is for mixed flow to occur where deformation 
rates have components parallel and perpendicular to the principal flow direction. 
In such flows, the elongational components may be associated with the converging- 
diverging flow paths [3, 24]. 

A general consensus has emerged recently that the flow through packed beds has 
a substantial extensional component and typical polymer solutions exhibit strain 
hardening in extension, which is mainly responsible for the reported dramatic in- 
creases in pressure drop. Thus in principle the shear viscosity alone is inadequate to 
explain the observed excessive pressure gradients. It is therefore essential to know 
the relative importance of elastic and viscous effects or equivalently the relationship 
between normal and shear stresses for different shear rates [10, 15, 25]. 

Elongational flow is fundamentally different from shear, the material property 
characterizing the flow is not the viscosity, but the elongational viscosity. The 
behavior of the extensional viscosity function is very often qualitatively different 
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from that of the shear viscosity function. For example, highly elastic polymer 
solutions that possess a viscosity which decreases monotonically in shear often 
exhibit an extensional viscosity that increases dramatically with strain rate. Thus, 
while the shear viscosity is shear-thinning the extensional viscosity is extension 
thickening [3, 5, 13]. 

3.2 Converging-Diverging Geometry 

An important aspect that characterizes the flow in porous media and makes it 
distinct from bulk is the presence of converging- diverging flow paths. This geo- 
metric factor significantly affects the flow and accentuate elastic responses. Several 
arguments have been presented in the literature to explain the effect of converging- 
diverging geometry on the flow behavior. Despite this diversity, there is a general 
consensus that in porous media the converging-diverging nature of the flow paths 
brings out both the extensional and shear properties of the fluid. The principal 
mode of deformation to which a material element is subjected as the flow converges 
into a constriction involves both a shearing of the material element and a stretch- 
ing or elongation in the direction of flow, while in the diverging portion the flow 
involves both shearing and compression. The actual channel geometry determines 
the ratio of shearing to extensional contributions. In many realistic situations in- 
volving viscoelastic flows the extensional contribution is the most important of 
the two modes. As porous media flow involves elongational flow components, the 
coil-stretch phenomenon can also take place. Consequently, a suitable model pore 
geometry is one having converging and diverging sections which can reproduce the 
elongational nature of the flow, a feature that is not exhibited by straight capillary 
tubes [10, 16, 23, 26, 27]. 

For long time, the straight capillary tube has been the conventional model for 
porous media and packed beds. Despite the general success of this model with 
the Newtonian and inelastic non-Newtonian flow, its failure with elastic flow was 
remarkable. To redress the flaws of this model, the undulating tube and converging- 
diverging channel were proposed in order to include the elastic character of the 
flow. Various corrugated tube models with different simple geometries have been 
used as a first step to model the effect of converging-diverging geometry on the 
flow of viscoelastic fluids in porous media (e.g. [16, 22, 28]). Those geometries 
include conically shaped sections, sinusoidal corrugation and abrupt expansions 
and contractions. Similarly, a bundle of converging-diverging tubes forms a better 
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model for a porous medium in viscoelastic flow than the normally used bundle of 
straight capillary tubes, as the presence of diameter variations makes it possible to 
account for elongational contributions. 

Many investigators have attempted to capture the role of the successive converging- 
diverging character of packed bed flow by numerically solving the flow equations 
in conduits of periodically varying cross sections. Different opinions on the success 
of these models can be found in the literature. Examples include [10, 15, 23, 26, 
29, 30]. With regards to modeling viscoelastic flow in regular or random networks 
of converging-diverging capillaries, very little work has been done. 
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4 Modeling the Flow in Porous Media 

The basic equations describing the flow of fluids consist of the basic laws of con- 
tinuum mechanics which are the conservation principles of mass, energy and linear 
and angular momentum. These governing equations indicate how the mass, energy 
and momentum of the fluid change with position and time. The basic equations 
have to be supplemented by a suitable rheological equation of state, or constitutive 
equation describing a particular fluid, which is a differential or integral mathemat- 
ical relationship that relates the extra stress tensor to the rate-of-strain tensor in 
general flow condition and closes the set of governing equations. One then solves 
the constitutive model together with the conservation laws using a suitable method 
to predict velocity and stress fields of the flows [7, 12, 25, 31, 32, 33]. 

The mathematical description of the flow in porous media is extremely com- 
plex task and involves many approximations. So far, no analytical fluid mechanics 
solution to the flow through porous media has been found. Furthermore, such a 
solution is apparently out of reach for the foreseeable future. Therefore, to investi- 
gate the flow through porous media other methodologies have been developed, the 
main ones are the macroscopic continuum approach and the pore-scale numerical 
approach. 

The widely used continuum models represent a simplified macroscopic approach 
in which the porous medium is treated as a continuum. All the complexities and 
fine details of the microscopic pore structure are absorbed into bulk terms like 
permeability that reflect average properties of the medium. Semi-empirical equa- 
tions such as Darcy's law, Blake-Kozeny-Carman or Ergun equation fall into this 
category. Several continuum models are based in their derivation on the capillary 
bundle concept. The advantage of the continuum method is that it is simple and 
easy to implement with no computational cost. The disadvantage is that it does not 
account for the detailed physics at the pore level. One consequence of this is that 
in most cases it can only deal with steady-state situations with no time- dependent 
transient effects. 

In the numerical approach, a detailed description of the porous medium at pore- 
scale level is adopted and the relevant physics of flow at this level is applied. To find 
the solution, numerical methods, such as finite volume and finite difference, usually 
in conjunction with computational implementation are used. The advantage of the 
numerical method is that it is the most direct approach to describe the physical 
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situation and the closest to full analytical solution. It is also capable, in principle 
at least, to deal with time-dependent transient situations. The disadvantage is that 
it requires a detailed pore space description. Moreover, it is usually very complex 
and hard to implement and has a huge computational cost with serious convergence 
difficulties. Due to these complexities, the flow processes and porous media that 
are currently within the reach of numerical investigation are the most simple ones. 

Pore-scale network modeling is a relatively novel method developed to deal 
with flow through porous media. It can be seen as a compromise between these 
two extreme approaches as it partly accounts for the physics and void space de- 
scription at the pore level with reasonable and generally affordable computational 
cost. Network modeling can be used to describe a wide range of properties from 
capillary pressure characteristics to interfacial area and mass transfer coefficients. 
The void space is described as a network of pores connected by throats. The pores 
and throats are assigned some idealized geometry, and rules which determine the 
transport properties in these elements are incorporated in the network to compute 
effective transport properties on a mesoscopic scale. The appropriate pore-scale 
physics combined with a geologically representative description of the pore space 
gives models that can successfully predict average behavior [34, 35]. 

In our investigation to the flow of Bautista-Manero fluids in porous media we 
use network modeling. Our model uses three-dimensional networks built from a 
topologically-equivalent three-dimensional voxel image of the pore space with the 
pore sizes, shapes and connectivity reflecting the real medium. Pores and throats 
are modeled as having triangular, square or circular cross-section by assigning a 
shape factor which is the ratio of the area to the perimeter squared and obtained 
from the pore space image. Most of the network elements are not circular. To ac- 
count for the non-circularity when calculating the volumetric flow rate analytically 
or numerically for a cylindrical capillary, an equivalent radius Reg is defined: 



where the geometric conductance, G, is obtained empirically from numerical simu- 
lation. Two networks obtained from Statoil and representing two different porous 
media have been used: a sand pack and a Berea sandstone. These networks are 
constructed by 0ren and coworkers [36, 37] from voxel images generated by sim- 
ulating the geological processes by which the porous medium was formed. The 
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physical and statistical properties of the networks are given in [38]. 

Assuming a laminar, isothermal and incompressible flow at low Reynolds num- 
ber, the only equations that need to be considered are the constitutive equation 
for the particular fluid and the conservation of volume as an expression for the 
conservation of mass. Because initially the pressure drop in each network element 
is not known, an iterative method is used. This starts by assigning an effective vis- 
cosity to each network element. The effective viscosity is defined as that viscosity 
which makes PoiseuUe's equation fit any set of laminar flow conditions for time- 
independent fluids [1]. By invoking the conservation of volume for incompressible 
fluid, the pressure field across the entire network is solved using a numerical solver 
[39]. Knowing the pressure drops, the effective viscosity of each element is updated 
using the expression for the flow rate with a pseudo-Poiseuille definition. The 
pressure field is then recomputed using the updated viscosities and the iteration 
continues until convergence is achieved when a specified error tolerance in total 
flow rate between two consecutive iteration cycles is reached. Finally, the total 
volumetric flow rate and the apparent viscosity in porous media, defined as the 
viscosity calculated from the Darcy's law, are obtained. 

With regards to modeling the flow in porous media of complex fluids which 
have time dependency due to thixotropic or elastic nature, there are three major 
difficulties : 

• The difficulty of tracking the fluid elements in the pores and throats and 
identifying their deformation history, as the path followed by these elements 
is random and can have several unpredictable outcomes. 

• The mixing of fluid elements with various deformation history in the individ- 
ual pores and throats. As a result, the viscosity is not a well-defined property 
of the fluid in the pores and throats. 

• The change of viscosity along the streamline since the deformation history is 
constantly changing over the path of each fluid element. 

In the current work, we deal only with one case of steady-state viscoelastic flow, 
and hence we did not consider these complications in depth. Consequently, the 
tracking of fluid elements or flow history in the network and other dynamic as- 
pects are not implemented in the non-Newtonian code. However, time-dependent 
effects in steady-state conditions are accounted for in the Tardy algorithm which 
is implemented in the code and will be presented in detail in § (5.1). 
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In our modeling approach, to solve the pressure field across a network of n 
nodes we write n equations in n unknowns which are the pressure values at the 
nodes. The essence of these equations is the continuity of fiow of incompressible 
fiuid at each node in the absence of source and sink. We solve this set of equations 
subject to the boundary conditions which are the pressures at the inlet and outlet. 
This unique solution is 'consistent' and 'stable' as it is the only mathematically 
acceptable solution to the problem, and, assuming the modeling process and the 
mathematical technicalities are correct, should mimic the unique physical reality 
of the pressure field in the porous medium. 
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5 Bautista-Manero Model 



This is a relatively simple model that combines the Oldroyd-B constitutive equation 
for viscoelasticity and the Fredrickson's kinetic equation for flow-induced structural 
changes usually associated with thixotropy. The model requires six parameters that 
have physical significance and can be estimated from rheological measurements. 
These parameters are the low and high shear rate viscosities, the elastic modulus, 
the relaxation time, and two other constants describing the build up and break 
down of viscosity. 

The Oldroyd-B model is a simplification of the more elaborate and rarely used 
Oldroyd 8-constant model which also contains the upper convected, the lower con- 
vected, and the corotational Maxwell equations as special cases. Oldroyd-B is the 
second simplest nonlinear viscoelastic model and is apparently the most popular 
in viscoelastic flow modeling and simulation. It is the nonlinear equivalent of the 
linear Jeffreys model, and hence it takes account of frame invariance in the non- 
linear regime. Consequently, in the linear viscoelastic regime the Oldroyd-B model 
reduces to the linear Jeffreys model. The Oldroyd-B model can be obtained by 
replacing the partial time derivatives in the differential form of the Jeffreys model 
with the upper convected time derivatives [12] 



where r is the extra stress tensor, Ai is the relaxation time, A2 is the retardation 
time, fio is the low-shear viscosity, 7 is the rate-of-strain tensor, and r is the upper 
convected time derivative of the stress tensor: 



where t is time, v = {vx,Vy,Vz) is the fluid velocity vector, (■) is the transpose of 
the tensor and Vv is the fluid velocity gradient tensor defined by 




(2) 




(3) 



Vv 



dvx dvx dvx 

dx dy dz 

dVy dVy dVy 

dx dy dz 



(4) 



, dvz dvz dvz J 
\ dx dy dz / 



V 

Similarly, 7 is the upper convected time derivative of the rate-of-strain tensor 
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given by: 



V d'j 



+ V ■ V7 — (Vv) ■ 7 — 7 ■ Vv 



(5) 




The kinetic equation of Fredrickson that accounts for the destruction and con- 
struction of structure is given by 



where /i is the non-Newtonian viscosity, t is the time of deformation, A is the 
relaxation time upon the cessation of steady flow, fio and fioo are the viscosities 
at zero and infinite shear rates respectively. A; is a parameter that is related to a 
critical stress value below which the material exhibits primary creep, r is the stress 
tensor and 7 is the rate of strain tensor. In this model, A is a structural relaxation 
time, whereas k is a kinetic constant for structure break down. The elastic modulus 
Go is related to these parameters by Go = [40, 41, 42, 43]. 

Bautista-Manero model was originally proposed for the rheology of worm-like 
micellar solutions which usually have an upper Newtonian plateau, and show strong 
signs of shear-thinning. The model, which incorporates shear-thinning, elasticity 
and thixotropy, can be used to describe the complex rheological behavior of vis- 
coelastic systems that also exhibit thixotropy and rheopexy under shear flow. The 
model predicts creep behavior, stress relaxation and the presence of thixotropic 
loops when the sample is subjected to transient stress cycles. The Bautista-Manero 
model has also been found to fit steady shear, oscillatory and transient measure- 
ments of viscoelastic solutions [40, 41, 42, 43]. 

5.1 Tardy Algorithm 

This algorithm is proposed by Philippe Tardy to compute the pressure drop-flow 
rate relationship for the steady-state flow of a Bautista-Manero fluid in simple 
capillary network models. The bulk rheology of the fluid and the dimensions of 
the capillaries making up the network are used as inputs to the models [43] . In the 
following paragraphs we outline the basic components of this algorithm and the 
logic behind it. This will be followed by some mathematical and technical details 
related to the implementation of this algorithm in our non-Newtonian code. 
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The flow in a single capillary can be described by the following general relation 

Q = CAP (7) 

where Q is the volumetric flow rate, G' is the flow conductance and AP is the 
pressure drop. For a particular capillary and a specific fluid, G' is given by 

G' = G'{fi) = constant Newtonian Fluid 

G' = G'{fi,AP) Purely viscous non-Newtonian Fluid 

G' = G'{fi,AP,t) Fluid with memory (8) 

For a network of capillaries, a set of equations representing the capillaries and 
satisfying mass conservation should be solved simultaneously to produce a consis- 
tent pressure field, as presented in § (4). For Newtonian fluid, a single iteration 
is needed to solve the pressure field since the conductance is known in advance as 
the viscosity is constant. For purely viscous non-Newtonian fluid, we start with an 
initial guess for the viscosity, as it is unknown and pressure-dependent, and solve 
the pressure field iteratively updating the viscosity after each iteration cycle until 
convergence is reached. For memory fluids, the dependence on time must be taken 
into account when solving the pressure field iteratively. Apparently, there is no 
general strategy to deal with such situation. However, for the steady-state flow of 
memory fluids a sensible approach is to start with an initial guess for the flow rate 
and iterate, considering the effect of the local pressure and viscosity variation due 
to converging-diverging geometry, until convergence is achieved. This approach is 
adopted by Tardy to find the flow of a Bautista-Manero fluid in a simple capil- 
lary network model. In general terms, the Tardy algorithm can be summarized as 
follows [43] 

• For fluids without memory the capillary is unambiguously defined by its ra- 
dius and length. For fluids with memory, where going from one section to 
another with different radius is important, the capillary should be modeled 
with contraction to account for the effect of converging-diverging geometry 
on the flow. The reason is that the effects of fluid memory take place on 
going through a radius change, as this change induces a change in shear 
rate and generation of extensional flow fields with a consequent viscoelastic 
and thixotropic effects. Examples of the converging- diverging geometries are 
given in Figure (17). 
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• Each capillary is discretized in the flow direction and a discretized form of the 
flow equations is used assuming a prior knowledge of the stress and viscosity 
at the inlet of the network. 

• Starting with an initial guess for the flow rate and using iterative technique, 
the pressure drop as a function of the flow rate is found for each capillary. 

• Finally, the pressure field for the whole network is found itcratively until 
convergence is achieved. Once the pressure field is found the flow rate through 
each capillary in the network can be computed and the total flow rate through 
the network can be determined by summing and averaging the flow through 
the inlet and outlet capillaries. 

A modified version of the Tardy algorithm was implemented in our non-Newtonian 
code. In the following, we outline the main steps of this algorithm and its imple- 
mentation. 



• From a one-dimensional steady-state Fredrickson equation in which the par 

d d 

dt dx 



tial time derivative is written in the form ^ = V-k-., the following equation 



can be obtained 



V^^^(^^\+kJ'^^]r^ (9) 



In this formulation, the fluid speed V is given by Q/nr'^ and the average shear 
rate in the tube with radius r is given by Q/nr^. 

• By similar argument, another simplified equation can be obtained from the 
Oldroyd-B model 

• The converging-diverging feature of the capillary is implemented in the form 
of a parabolic profile as outlined in Appendix A. 



• Each capillary in the network is discretized into m slices, each with width 
5x — L/m where L is the capillary length. 
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From Equation (10), applying the simplified assumption ^ = ^'^^Jj^^ where 
the subscripts 1 and 2 stand for the inlet and outlet of the slice respectively, 
we obtain an expression for T2 in terms of Ti and fi2 



2/i27 + 



^2 = -Wv^ (11) 



1 + 



GoSx 



From Equations (9) and (11), applying ^ = ^'^^^^^^ , a third order polynomial 
in /i2 is obtained. The coefficients of the four terms of this polynomial are 



V 2kY k^nV 



XfloGoSx floe flooGoSx 



V 1 



1 -U i -I- ^ ^1 

^2 ■ A Go{6xy 

A: 

The algorithm starts by assuming a Newtonian flow in a network of straight 
capillaries. Accordingly, the volumetric flow rate Q, and consequently V and 
7, for each capillary are obtained. 

Starting from the inlet of the network where the viscosity and the stress are 
assumed to have known values of /io and RAP/2L for each capillary respec- 
tively, the computing of the non-Newtonian flow in a converging-diverging 
geometry takes place in each capillary independently by calculating and 
T2 slice by slice, where the values from the previous shce are used for /zi and 
Ti of the current slice. 

For the capillaries which are not at the inlet of the network, the initial values 
of the viscosity and stress at the inlet of the capillary are found by computing 
the Q-weighted average from all the capillaries that feed into the pore which 
is at the inlet of the corresponding capillary. 

To flnd /i2 of a slice, a bisection numerical method is used. To eliminate 
possible non-physical roots, the interval for the accepted root is set between 
zero and 3/io with /io used in the case of failure. These conditions are logical 
as long as the slice is reasonably thin and the flow and fluid are physically 
viable. In the case of a convergence failure, error messages are issued to 
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inform the user. No failure has been detected during the many runs of this 
algorithm. Moreover, extensive sample inspection of the /i2 values has been 
carried out and proved to be sensible and realistic. 

• The value found for the /i2 is used in conjunction with Equation (11) to find 
T2 which is needed as an input to the next slice. 

• Averaging the value of /xi and fi2, the viscosity for the slice is found and used 
with Poiseuille law to find the pressure drop across the slice. 

• The total pressure drop across the whole capillary is computed by summing 
up the pressure drops across its individual slices. This total is used with 
Poiseuille law to find the effective viscosity for the capillary as a whole. 

• Knowing the effective viscosities for all capillaries of the network, the pressure 
field is solved iteratively using an algebraic multi-grid solver, and hence the 
total volumetric fiow rate from the network and the apparent viscosity are 
found. 

5.2 Initial Results of the Modified Tardy Algorithm 

The modified Tardy algorithm was tested and assessed. Various qualitative aspects 
were verified and proved to be correct. Some general results and conclusions are 
outlined below with sample graphs and data for the sand pack and Berea net- 
works using a calculation box with x =0.5 and x =0.95. We would like to remark 
that despite our effort to use typical values for the parameters and variables, in 
some cases we were forced, for demonstration purposes, to use eccentric values to 
accentuate the features of interest. In Table (1) we present some values for the 
Bautista-Manero model parameters as obtained from the wormlike micellar system 
studied by Anderson and co-workers [44] to give an idea of the parameter ranges 
in real fluids. The system is a solution of a surfactant concentrate [a mixture of 
the cationic surfactant erucyl bis(hydroxyethyl)methylammonium chloride (EHAC) 
and 2-propanol in a 3:2 weight ratio] in an aqueous solution of potassium chloride. 

5.2.1 Convergence-Divergence 

A dilatant effect due to the converging- diverging feature has been detected relative 
to a network of straight capillaries. The viscosity increase is a natural viscoelastic 
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Table 1: Some values of the wormlike micellar system studied by Anderson and 
co-workers [44] which is a solution of surfactant concentrate (a mixture of EHAC 
and 2-propanol) in an aqueous solution of potassium chloride. 



Parameter 


Value 


Go (Pa) 


1 - 10 


fio (Pa.s) 


115 - 125 


/ioo (Pa.s) 


0.00125 - 0.00135 


A(s) 


1 - 30 


k (Pa-i) 


10"' - 10"' 



response to the radius tightening at the middle. As the corrugation feature of 
the tubes is exacerbated by narrowing the radius at the middle, the viscoelastic 
dilatant behavior is intensified. The natural explanation is that the increase in 
apparent viscosity should be proportionate to the magnitude of tightening. This 
feature is presented in Figures (1) and (2) for the sand pack and Berea networks 
respectively on a log-log scale. 

5.2.2 Diverging-Converging 

The investigation of the effect of diverging- converging geometry by expanding the 
radius of the capillaries at the middle revealed a thinning effect relative to the 
straight and converging-diverging geometries, as seen in Figures (3) and (4) for 
the sand pack and Berea networks respectively on a log-log scale. This is due 
to viscoelastic response in the opposite sense to that of converging-diverging by 
enlarging the radius at the middle. 

5.2.3 Number of Slices 

As the number of slices of the capillaries increases, the algorithm converges to a 
stable and constant solution within acceptable numerical errors. This indicates 
that the numerical aspects of the algorithm are functioning correctly because the 
effect of discretization errors is expected to diminish by increasing the number 
of slices and hence decreasing their width. A sample graph of apparent viscosity 
versus number of slices for a typical data point is presented in Figure (5) for the 
sand pack and in Figure (6) for Berea. 



5.2.3 Number of Slices 



24 




1.0E-06 1.0E-05 1.0E-04 1.0E-03 1.0E-02 1.0E-01 

Darcy Velocity / (m/s) 

Figure 1: The Tardy algorithm sand pack results for Go=0.1Pa, //oo=0.001 Pa.s, 
/io=1.0Pa.s, A=1.0s, k^lO~^ Pa.~^, /e=1.0, m=10 slices, with varying (1.0, 0.8, 
0.6 and 0.4). 

1.0E+01 T 




1.0E-06 1.0E-05 1.0E-04 1.0E-03 1.0E-02 

Darcy Velocity / (m/s) 

Figure 2: The Tardy algorithm Bcrea results for Go=0.1Pa, /ioo=0.001 Pa.s, 
/Xo=1.0Pa.s, A=1.0s, A;=10~^Pa~^, /e=1.0, m=10 slices, with varying (1.0, 
0.8, 0.6 and 0.4). 
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1.0E-05 1.0E-04 1.0E-03 1.0E-02 1.0E-01 

Darcy Velocity / (m/s) 

Figure 3: The Tardy algorithm sand pack results for Go=0.1Pa, //oo=0.001 Pa.s, 
/io=l-0 Pa.s, A=1.0s, k—10~^ Pa~^, fe—1.0, m—10 slices, with varying fj^ (0.8, 1.0, 
1.2, 1.4 and 1.6). 
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1.0E-06 1.0E-05 1.0E-04 1.0E-03 1.0E-02 
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Figure 4: The Tardy algorithm Berea results for Go=0.1Pa, /ioo=0. 001 Pa.s, 
/io=1.0Pa.s, A=1.0s, k—10~^Pa,~^, fe—1.0, m=10 slices, with varying (0.8, 
1.0, 1.2, 1.4 and 1.6). 

5.2.4 Boger Fluid 



A Boger fluid behavior was observed when setting fio = /^oo- However, the apparent 
viscosity increased as the converging-diverging feature is intensified. This feature is 
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Figure 5: The Tardy algorithm sand pack results for Pa, yUoo=0.001 Pa.s, 

/io=l-OPa.s, A=1.0s, fc=10~^Pa~^, /e=1.0, /m=0.5, with varying number of slices 
for a typical data point (AP=100Pa). 
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Figure 6: The Tardy algorithm Berea results for Go=l-OPa, /ioo=0.001 Pa.s, 
yUo=l-OPa.s, A=1.0s, /c=10~^Pa~^, /e=1.0, /m=0.5, with varying number of slices 
for a typical data point (AP=200Pa). 

demonstrated in Figure (7) for the sand pack and in Figure (8) for Berea on a log- 
log scale. Since Boger fluid is a limiting and obvious case, this behavior indicates 



5.2.5 Elastic Modulus 



27 



that the model, as implemented, is well-behaved. The viscosity increase is a natural 
viscoelastic response to the radius tightening at the middle, as discussed earlier. 

5.2.5 Elastic Modulus 

The effect of the elastic modulus Go was investigated for shear-thinning fluids, 
i.e. /io > fioo, by varying this parameter over several orders of magnitude while 
holding the others constant. It was observed that by increasing Go, the high- 
shear apparent viscosities were increased while the low-shear viscosities remained 
constant and have not been affected. However, the increase at high-shear rates has 
almost reached a saturation point where beyond some limit the apparent viscosities 
converged to certain values despite a large increase in Go- A sample of the results for 
this investigation is presented in Figure (9) for the sand pack and in Figure (10) for 
Berea on a log-log scale. The stability at low-shear rates is to be expected because 
at low flow rate regimes near the lower Newtonian plateau the non-Newtonian 
effects due to elastic modulus are negligible. As the elastic modulus increases, an 
increase in apparent viscosity due to elastic effects contributed by elastic modulus 
occurs. Eventually, a saturation will be reached when the contribution of this factor 
is controlled by other dominant factors and mechanisms from the porous medium 
and flow regime. 

5.2.6 Shear-Thickening 

A slight shear-thickening effect was observed when setting fio < /ioo while holding 
the other parameters constant. The effect of increasing Go on the apparent vis- 
cosities at high-shear rates was similar to the effect observed for the shear-thinning 
fluids though it was at a smaller scale. However, the low-shear viscosities are not 
affected, as in the case of shear-thinning fluids. A convergence for the appar- 
ent viscosities at high-shear rates for large values of Go was also observed as for 
shear-thinning. A sample of the results obtained in this investigation is presented 
in Figures (11) and (12) for the sand pack and Berea networks respectively on a 
log-log scale. It should be remarked that as the Bautista-Manero is originally a 
shear-thinning model, it should not be expected to fuUy-predict shear-thickening 
phenomenon. However, the observed behavior is a good sign for our model be- 
cause as we push the algorithm beyond its limit, the results are still qualitatively 
reasonable. 
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Figure 7: The Tardy algorithm sand pack results for Go=1.0Pa, //oo=A*o=0.1 Pa.s, 
A=1.0s, A;=10-5Pa-\ /e=1.0, m=10 slices, with varying (0.6, 0.8, 1.0, 1.2 and 
1.4). 
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Figure 8: The Tardy algorithm Berea results for 6*0=1-0 Pa, /ioo=/^o=0.1 Pa.s, 
A=1.0s, k=10-^Pci-\ /e=1.0, m=10 slices, with varying (0.6, 0.8, 1.0, 1.2 
and 1.4). 
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Figure 9: The Tardy algorithm sand pack results for //oo=0.001 Pa.s, /Xo=l-OPa.s, 

A=1.0s, A':=10-5Pa-^ /e=1.0, m=10 slices, with varying (0.1, 1.0, 10 
and 100 Pa). 
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Figure 10: The Tardy algorithm Berea results for /Xoo=0.001 Pa.s, /io= l.OPa.s, 
A=1.0s, A;=10-5Pa-^ /e=1.0, /^=0.5, m=10 slices, with varying Go (0.1, 1.0, 10 
and 100 Pa). 

5.2.7 Relcixation Time 



The effect of the structural relaxation time A was investigated for shear-thinning 
fluids by varying this parameter over several orders of magnitude while holding 
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Figure 11: The Tardy algorithm sand pack results for iJ,^—10.0Pa,.s, /Xo=0.1Pa.s, 
A=1.0s, /c=10~^Pa~^, /e=1.0, fjn—0.5, m=10 slices, with varying Go (0.1 and 
100 Pa). 




Figure 12: The Tardy algorithm Berea results for /ioo=10.0 Pa.s, /io=0.1Pa.s, 
A=1.0s, A;=10~*^ Pa~^, /e=1.0, /m=0.5, m=10 slices, with varying Go (0.1 and 
100 Pa). 



the others constant. It was observed that by increasing the structural relaxation 
time, the apparent viscosities were steadily decreased. However, the decrease at 
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high-shear rates has almost reached a saturation point where beyond some hmit 
the apparent viscosities converged to certain values despite a large increase in A. 
A sample of the results is given in Figure (13) for the sand pack and Figure (14) 
for Berea on a log-log scale. The effects of relaxation time are expected to ease 
as the relaxation time increases beyond a limit such that the effect of interaction 
with capillary constriction is negligible. Despite the fact that this feature requires 
extensive investigation for quantitative confirmation, the observed behavior seems 
qualitatively reasonable considering the flow regimes and the size of relaxation 
times of the sample results. 

5.2.8 Kinetic Parameter 

The effect of the kinetic parameter for structure break down k was investigated for 
shear-thinning fluids by varying this parameter over several orders of magnitude 
while holding the others constant. It was observed that by increasing the kinetic 
parameter, the apparent viscosities generally decreased. However, in some cases 
the low-shear viscosities has not been substantially affected. A sample of the results 
is given in Figures (15) and (16) for the sand pack and Berea networks respectively 
on a log-log scale. The decrease in apparent viscosity on increasing the kinetic 
parameter is a natural response as the parameter quantifies structure break down 
and hence reflects thinning mechanisms. It is a general trend that the low flow 
rate regimes near the lower Newtonian plateau are usually less influenced by non- 
Newtonian effects. It will therefore be natural that the low-shear viscosities in 
some cases experienced minor changes. 
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Figure 13: The Tardy algorithm sand pack results for Go=1.0Pa, //oo=0.001 Pa.s, 
;Uo=1.0Pa.s, A;=10~^Pa~^, fe—l.O, /^^O.S, m—10 slices, with varying A (0.1, 1.0, 
10, and 100 s). 
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Figure 14: The Tardy algorithm Berea results for Go=1.0Pa, /ioo=0.001 Pa.s, 
/Xo=l-OPa.s, A;=10~^Pa~\ /e=l-0, /^=0.5, m=10 slices, with varying A (0.1, 1.0, 
10, and 100 s). 
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Figure 15: The Tardy algorithm sand pack results for Go=1.0Pa, //oo=0.001 Pa.s, 
//o=l-OPa.s, A=10s, /e=1.0, m—10 slices, with varying k (10~^, 10~^, 

10-^ 10-6 and 10"^ Pa"^). 
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Figure 16: The Tardy algorithm Berea results for Go=1.0Pa, /ioo=0.001 Pa.s, 
/Xo=l-OPa.s, A=1.0s, fe—1.0, fjn—0.5, m—10 slices, with varying k (10~^, 10~^, 
10-^ 10-6 and 10"^ Pa-^). 



6 Conclusions 



A non-Newtonian flow simulation code based on pore-scale network modeling has 
been extended to account for viscoelastic and thixotropic behavior using a Bautista- 
Manero fluid model. The basis of the implementation is a numerical algorithm 
originally suggested by Philippe Tardy. A modified version of this algorithm was 
used to investigate the effect of converging-diverging geometry on the steady- 
state flow. The implementation of this model has been examined and assessed 
using a sand pack and a Berea networks. The generic behavior indicates qual- 
itatively correct trends. A conclusion was reached that the current implemen- 
tation is a proper basis for investigating the steady-state viscoelastic effects and 
some thixotropic effects due to converging-diverging geometries in porous media. 
This implementation can be used as a suitable foundation for further develop- 
ment. The future work can include elaborating the modified Tardy algorithm and 
thoroughly investigating its viscoelastic and thixotropic predictions in quantitative 
terms. The code (called Non-Newtonian Code 2) with complete documentation 
and the networks can be obtained from this URL: http : //www3 . imperial . ac . uk/ 
earthsciencecindengineerinresearch/perm/porescalemodellinsof tware/non-newtoiiiany„ 
20code. 
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Nomenclature 



7 strain rate (s^^) 

7 ratc-of-strain tensor 

A structural relaxation time in Fredrickson model (s) 

Ai relaxation time (s) 

A2 retardation time (s) 

II viscosity (Pa.s) 

/lo zero-shear viscosity (Pa.s) 

//oo infinite-shear viscosity (Pa.s) 

T stress (Pa) 

T stress tensor 

fe scale factor for the entry of corrugated tube ( — ) 

f„i scale factor for the middle of corrugated tube ( — ) 

G geometric conductance (m^) 

G' flow conductance (m'^.Pa~^.s~^) 

Go elastic modulus (Pa) 

k parameter in Fredrickson model (Pa~^) 

L tube length (m) 

P pressure (Pa) 

AP pressure drop (Pa) 

Q volumetric flow rate (m^.s~^) 

r radius (m) 

R tube radius (m) 

Req equivalent radius (m) 

Rmax maximum radius of corrugated capillary 

Rmin minimum radius of corrugated capillary 

t time (s) 

V fluid velocity vector 

Vv fluid velocity gradient tensor 
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V fluid speed (m.s ^) 

5x smaU change in x (m) 



upper convected time derivative 
(•)'^ matrix transpose 

Xi network lower boundary in the non-Newtonian code 

x^ network upper boundary in the non-Newtonian code 

Note: units, when relevant, are given in the SI system. Vectors and tensors are 
marked with boldface. Some symbols may rely on the context for unambiguous 
identification. 
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Appendix A: Converging-Diverging Geom- 
etry and Tube Discretization 



There are many simplified converging-diverging geometries for corrugated capillar- 
ies which can be used to model the flow of viscoelastic fluids in porous media; some 
suggestions are presented in Figure (17). In this study we adopted the paraboloid 
and hence developed simple formulae to find the radius as a function of the distance 
along the tube axis, assuming cylindrical coordinate system, as shown in Figure 
(18). 

For the paraboloid depicted in Figure (18), we have 

r{z) = az^ + bz + c (13) 

Since 



r{-L/2)=r{L/2) = Rma. & r(0) = (14) 

the paraboloid is uniquely identified by these three points. On substituting and 
solving these equations simultaneously, we obtain 



2^^ 
L 



that is 



{Rmax Rmin) b — & C — Rmin iX^) 



In the modified Tardy algorithm for steady-state viscoelastic flow as imple- 
mented in our non-Newtonian code, when a capillary is discretized into m slices 
the radius r[z) is sampled at m 2;-points given by 

z = ~ + k- (A; = l,...,m) (17) 

2 m 

More complex polynomial and sinusoidal converging-diverging profiles can also 
be modeled using this approach. 
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Figure 17: Examples of corrugated capillaries which can be used to model 
converging-diverging geometry in porous media. 



r 




Figure 18: Radius variation in the axial direction for a corrugated paraboloid 
capillary using a cylindrical coordinate system. 
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